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Simulation of a Langevin-dynamics model demonstrates emergence of critical fluctuations and 
anomalous grain transport which have been observed in experiments on "soft" quasi-two-dimensional 
dusty plasma clusters. Our model does not contain external drive or plasma interactions that serve 
to drive the system away from thermodynamic equilibrium. The grains are confined by an external 
potential, interact via static Yukawa forces, and are subject to stochastic heating and dissipation 
from neutrals. One remarkable feature is emergence of leptokurtic probability distributions of grain 
displacements £(t) on time-scales r < ta , where ta is the time at which the standard deviation 
<t(t) = {£, 2 (t)) 1 ^ 2 approaches the mean inter-grain distance A. Others are development of humps 
in the distributions on multiples of A, anomalous Hurst exponents, and transitions from leptokurtic 
towards Gaussian displacement distributions on time scales r > ta- The latter is a signature of 
intermittency, here interpreted as a transition from bursty transport associated with hopping on 
intermediate time scales to vortical flows on longer time scales. These intermittency features are 
quantitatively modeled by a single-partice Ito-Langevin stochastic equation with a non-linear drift 
term. 



I. INTRODUCTION 

Since a large aggregate of dust grains embedded in a 
plasma can attain virtually all known states of matter 
[l| , the ability of optical tracking of individual grain tra- 
jectories has made these systems particularly attractive 
for studying transport properties and emergent complex 
behavior in soft materials. 

In a number of recent dusty plasmas experiments 
anomalous properties, like non-Gaussian displacement 
distributions and anomalous diffusion, have been ob- 
served. Yet it is not clear whether this behavior derives 
from non-equilibrium physics due to external forcing of 
the grain aggregate, or whether critical fluctuations Q 
and anomalous transport may emerge in aggregates sub- 
ject exclusively to statically screened Coulomb forces and 
independent stochastic forcing and friction due to colli- 
sions with the background neutral gas. The purpose of 
this paper is to demonstrate that the latter is possible, 
and that non-Gaussian statistics does not have to arise 
from particular non-equilibrium physics. 

So far few molecular-dynamics (MD) or Langevin- 
dynamics (LD) simulations of dust transport have been 
able to reproduce the anomalous results observed in ex- 
periments, although subdiffusion on short time scales 
have recently been reported for a bidispersive disordered 
assembley of Yukawa particles p|. Here we adopt the 
terminology employed in [J] , where MD refers to simula- 
tion of an aggregate of particles interacting via Yukawa 
forces and LD refers to a system that in addition is sub- 
ject to stochastic forcing and friction represented by a 
Langevin thermostat. The latter represents a heat bath 
which should allow the system to relax to thermodynamic 
equilibrium with the neutral gas. 

Under ground-based laboratory conditions dust grains 
levitate in electrostatic sheaths, which give rise to plasma 



anisotropy and supra-thermal ion flows that may act 
as external forcing. Such experiments are referred to 
as two-dimensional (2D or quasi-2D) due to the fact 
that vertical motion of particles is severely restricted 
(monolayers) and/or particles are vertically aligned due 
to the ion flow. Dust transport under such condi- 
tions has been extensively studied both in laboratory 
0, 0, II pjMMllMl- El El EE E3, El and simula- 
tions [J, [la R6 . [l7|, El I • Anomalous transport exhibi ting 
Gaussian ^.Tla] as well as non-Gaussian [E EE EE EE EE] 
statistics has been reported for both 2D and quasi-2D 
cases. 

In this paper we show that a model that includes only 
Yukawa interaction between particles, an externally im- 
posed confining potential, and stochastic heating and dis- 
sipation from neutrals, is able to demonstrate emergence 
of non-Gaussian statistical features associated with ob- 
served cooperative particle hopping 0, [E EH EE EH and 
vortical flow patterns [3, E3 • The simulations also re- 
veal some features for which we report the experimen- 
tal footprint here for the first time - the development 
of humps on the displacement probability distribution 
function (PDF) corresponding to integer multiples of the 
mean inter-grain distance A. 

We study quasi-2D systems with varying "stiffness", 
with particular attention paid to the transition in trans- 
port characteristics from time scales t <C ta to r > ta, 
where ta is the characteristic time for a dust grain to dif- 
fuse a distance A, i.e. ct(ta) = (^(ta)) 1 / 2 = A. For this 
purpose we analyze data from ground-based laboratory 
experiments (quasi-2D) and perform molecular dynamics 
simulations in quasi-2D and 2D over a range of states, 
from pure crystal to liquid with particular focus on soft 
(viscoelastic or hexatic) states at the transition between 
crystal and liquid. 

By examining and contrasting these data, the follow- 
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ing picture emerges: the shape of the position displace- 
ment distribution P[£(t)] is Gaussian for times shorter 
than the grain-grain collision time r c . On these short 
time scales the grain motion is ballistic, and the displace- 
ment divided by r can be interpreted as the instanta- 
neous grain velocity. Hence the Gaussian displacement 
PDF reflects the Maxwellian velocity distribution result- 
ing from the Gaussian stochastic heating term. On scales 
t c < t < Ta the PDF develops a leptokurtic shape. For 
stiffcr systems, as r approaches ta humps develop on the 
tail of the PDF at displacements £(t) corresponding to 
integer multiples of A. These humps signify that on these 
time scales the grain hopping in the imperfect hexagonal 
lattice influences the statistical transport characteristics. 
On even longer time scales, the kurtosis of the PDFs ap- 
proach that of a Gaussian, indicating that multiple hop- 
pings may represent a sequence of independent transport 
events @. 

In softer systems the picture on intermediate and 
longer time scales is different. The humps due to hopping 
are masked by the emergence of large-scale vortical mo- 
tions, and the PDFs on the hydrodynamical time-scales 
reflect these vortex motions. 

An interesting observation is that Gaussian statistics 
characterizes both the solidified crystalline state and the 
pure liquid state. Non-Gaussian statistics appears when 
the crystal is sufficiently soft to allow development of 
defects and particle hopping, but sufficiently stiff to re- 
strict free diffusion on scales larger that A. It appears on 
spatial displacement scales from a fraction of A and at 
least up to those scales where elasticity dominates over 
viscosity. The fact that Gaussian statistics and normal 
transport is restored in the purely liquid phase, indicates 
that the anomalous characteristics are associated with 
critical fluctuations in the form of a scale-free hierarchy 
of vortices The vortices emerge from the stochastic 
forcing only if the lattice resides in a soft, viscoelastic 
transitional state between crystal and liquid. 

We also discuss the possible physical origin of non- 
Gaussian displacement distributions and anomalous 
transport, and show how these feature derive naturally 
from a single-particle Langcvin equation with a non- 
linear drift (friction) term modeling the collective effects 
of dust-dust interaction. 



II. DATA SOURCES 
A. Laboratory experiments 

For this paper a re-analysis has been made of raw data 
from the experiment reported in Ref. fl3j . This experi- 
ment was performed in a capacitively coupled rf discharge 
operated in argon at pressure 1.8 Pa and rf power of 10 
W. A cluster of 120 melamine- formaldehyde monodis- 
persive spheres with a diameter of 8.9 /im, was confined 
by the potential created by a copper ring mounted on the 
lower electrode. The dust particles were illuminated by a 
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FIG. 1: (Color online) Grain trajectories from the molecular 
dynamics simulations of the Yukawa-Langevin system. The 
trajectories are color-coded by the time. Motion of particles 
during time 100 8t in the case of parabolic confinement (a) and 
the hard wall confinement (b). A typical trajectory from the 
outer region of panel (a) plotted for time 4000 St (c). Close- 
up on stiff central region of panel (b), motion of particles is 
shown for 100 St (d). 



horizontal laser sheet, and 50 000 images were taken by 
a video camera at a sampling rate of 30 Hz and spatial 
resolution of 6.7 /im/ pixel. 

As seen from Fig. 1(b) of Ref. [l3|, where trajectories 
of particles over time lag of 1.7 sec are presented, the 
grains have limited mobility, and we refer to these data as 
our "stiff cluster" experiment. The re-analysis shows the 
evolution of radial displacement PDFs computed with 
improved statistics compared to that of Ref. [l3| . This 
new analysis reveals features of the PDFs which were not 
resolved in the earlier data, and also allows us to obtain 
PDFs of displacements on much longer time scales. 

In Ref. [HI a large circular cluster of 600 dust grains in 
a monolayer configuration was studied in a similar exper- 
imental set-up. The cluster was brought to a partly or- 
dered state allowing large-scale viscoelastic vortex flows 
(see Figs. 1 and 2 of Ref. 14]). The experiment was 
operated in argon at a pressure 4 Pa and rf power of 19 
W with melamine-formaldchyde spheres with diameter of 
7.2 fim. The grains were confined radially by the poten- 
tial created by a cavity of 6 cm radius machined into the 
lower electrode. The cluster diameter was 16 mm with 
interparticle distance A sa 0.6 mm. In this experiment 
30 000 images were taken by a video camera at sampling 
rate of 30 Hz and spatial resolution of 24 /im/ 'pixel. We 
will refer to this data as the "soft cluster" experiment. 
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We will also connect the results from our model simula- 
tions to the experimental studies on cold quasi-2D dusty- 
plasma liquids reported in Refs. [7L ITTI. \l2\. 

B. Langevin dynamics simulations 

The simulations are performed for the Yukawa system 
of particles embedded into a Langevin thermostat [l9[ 
which simulates the effect of the collisions with the elec- 
trically neutral molecules of the background gas (details 
are given in [20I]). The pair interaction between the dust 
particles is described by the Debye-Hiickel (Yukawa) po- 
tential: 

<Kr) = (Q/r)exp(-r/A D ), (1) 

where Q = eZ,i is the particle charge and Ad is the 
the plasma screening length. Particle radius, mass and 
charge, mean inter-particle distance, neutral gas density 
and temperature are the parameters of the simulations. 
Initially N particles having the charge Q are randomly 
distributed over a prescribed simulation volume to reach 
a value of (average) screening parameter k = A/ Ad close 
to that in the experiments. After a transient phase, the 
system relaxes to a statistically stationary state with the 
dust velocity distribution in thermal equilibrium with the 
Langevin thermostat. This state does not depend on the 
initial conditions, which is also the case in the experi- 
ments. Three-dimensional particle dynamics of the i'th 
particle with mass m evolves according to the equation, 

TO ~y^ + mU "~d7 = ~® V ^ ~~ V<M»y) + Li(i) , 

(2) 

where v n = l/r„ is the neutral gas friction, <p c is the 
confinement potential and = | r j — Vj \ . The stochastic 
Langevin force Li(t) (thermal noise induced by neutral 
gas particles) is defined from the relation: 

(L,(t) ■ Lj(t + t))= 2u n mTS ij 6{r), (3) 

with the mean zero condition (Lj(t)) = 0, where T is the 
neutral gas temperature. 

Two types of the confinement have been used in the 
simulations; parabolic and hard wall. For the parabolic 
trap the confinement potential is 4*%( r ) K 7,2 ■ F° r the 
hard wall case the electric field associated with the con- 
finement exponentially increases to the boundaries, so 
that <fr% w (r) oc exp[(r — r b )/S w ] — 1 at r > r b and 
</>e W (r) = for r < r b , where rb is the size of the sys- 
tem and S w is a stiffness of the hard wall potential. In 
our simulations a value of S w ~ A/3 has been used. 

To simulate a quasi-2D configuration the motion in the 
vertical (z) direction has been suppressed by a strong 
confinement in this direction. In most of the runs per- 
formed, the particle displacements were restricted to a 
range smaller than the mean horizontal intcrparticle dis- 
tance A, which is the typical situation for monolayer ex- 
periments. In addition, simulations of 2D systems (only 



horizontal x- and y- components considered) have also 
been run for both types of the confinement. 

The parabolic confinement is in some respects a more 
realistic representation of the experiments discussed here 
and, as seen from Fig.[T](a) it reproduces flow fields quite 
similar to those observed in the experiments. As in the 
experiment, such confinement results in a nonuniform 
cluster with more packed and stiff core and looser bound- 
aries. The hard- wall potential, see Fig. QJb), yields a dif- 
ferent structure; particles are more packed at the bound- 
aries (for more details see 

HJIHIH). Ideally ' 

compar- 
ison of these different confinement schemes could help 
identifying transport features which are robust to varia- 
tions of the confinement. 

As one can observe from Fig. [3a) the resulting equi- 
librium state is a cluster of dust grains in a partially 
ordered (hexatic) phase, similar to those reported from 
2D Yukawa MD simulations on square double-periodic 
domains in e.g. [f| and in 2D hard-core simulations in 
Q. However, our cluster is of finite size (non-periodic 
boundary conditions), confined by a circular symmetric 
parabolic potential, and is radially nonuniform. Like in 
the experiments the simulations are intended to describe, 
the transport over large distances takes place through de- 
velopment and propagation of 5-fold and 7-fold topolog- 
ical crystal defects. 

By tuning the input parameters we obtain different 
states, ranging from pure crystal to pure liquid. As 
mentioned in the introduction we are mainly focused on 
semi-ordered (viscoclastic) states at the transition be- 
tween crystal and liquid. To obtain such states, the cou- 
pling parameter, defined as T s = e 2 Z%/(TA) cxp(— n), 
was changed in the range 50 — 100 by varying T . The 
latter also implies change of the inter-particle distance; 
the corresponding range of the screening parameter was 
k — 2, — 3, which is typical for the experiments. For each 
state studied, A was determined as the distance to the 
first local maximum of the pair correlation function av- 
eraged over the cluster. This distance varies less than 
20% from the core to the edge for parabolical confine- 
ment, and is nearly uniformly distributed in space for 
hard wall confinement. 

The particle positions have been sampled with a time 
resolution St chosen to be of the order of grain-grain col- 
lision time scale, St ~ r c . The grain-grain interaction 
time was estimated as t c ~ 1/(2itQe), where £Ie is the 
Einstein frequency which characterizes the linear oscilla- 
tions of individual particles in a lattice and is a function 
of the screening parameter n. For k ~ 2 — 3 the Einstein 
frequency is about a factor of two of the dust-plasma 
frequency Qq = \/e 2 Z 2 /mA 3 [l8[ . The dust-neutral col- 
lision time scale was r„ ~ 10<5t for all runs. 

For each state, a simulation of N = 1000 particles has 
been run during a time ~ 4000 St, thus proving adequate 
statistics for study of the system's evolution on long time 
scales. 

In both experiments and simulations the nature of the 
dynamics in the regimes where the system is in a semi- 
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FIG. 2: Experimental results from stiff cluster (Ref.JlJI). 
PDFs of radial displacements £ r (r) = r(t + r) — r(t) (dots) 
in units of interparticle distance, A, for different time lags; 
r = 0.033 s (a), r = 0.53 s (b), r = 8.5 s (c) and r = 136 s 
(d). The full curves are fits of the stretched Gaussian given 
by Eq. @ with q = 1.5 (a-d), and the thick full curve is a 
Gaussian fit (d). 



ordered state is very sensitive to fine tuning of the con- 
trol parameters, and we have not been able to reproduce 
all experimental features quantitatively. One reason for 
this may be that we are not able to fine tune parameters 
to the right state, but it can also be due to shortcom- 
ings of the model. However, the purpose of this paper 
is not to reproduce experimental results quantitatively, 
but to demonstrate that important qualitative features 
emerge from a simple Langevin model which does not in- 
volve other interactions between the dust particles and 
the plasma than the Debye screening of the Coulomb po- 
tential. 



III. RESULTS ON GRAIN TRANSPORT 

A. Laboratory results 

Time evolution of the PDFs of radial displacements 
£ r (r) = r(t + r) — r(t) in units of interparticle distance A 
for the stiff cluster is shown in Fig. [5] The distributions 
are presented on a logarithmic scale and the deviation 
from a Gaussian shape on small time scales Fig. [2] (a,b) 
is evident. On those scales the shape of the PDF is well 
described by the stretched Gaussian 



P(^ r ,r) = A(r)exp[-B(rMr\ p ] 



(4) 



where p = 1 corresponds to an exponential distribution 
and p = 2 to the Gaussian. For the data presented the 
nonlinear fit (thin full lines) yields the exponent p ss 1.5. 

A possible, and rather trivial, explanation of the 
stretched gaussian PDFs would be the spatial non- 
uniformity of the system. For instance, the dust grains 
in the laboratory clusters are more mobile close to the 
edge than in the central region. Thus, it is conceivable 



that a more heavy-tailed PDF could emerge from com- 
bining several Gaussian PDFs with different variances. 
We have tested this in the laboratory data by produc- 
ing separate PDFs for particle displacements in the core 
and edge regions, respectively. However, although these 
PDFs have different variances, they both retain essen- 
tially the same strctched-Gaussian PDF, and hence this 
is an intrinsic property of the particle dynamics on the 
short scales both in core an the edge, and is not an arti- 
fact due to spatial non-uniformity. 

When a significant number of particle displacements 
approach the mean inter-grain distance A the PDF starts 
to develop "shoulders" around A (Fig. Hfc)). The tran- 
sition towards a Gaussian (Fig. H{d), thick line) on the 
longest time scales involves a regime where the PDF de- 
velops "humps" at £ r nA, where n is a positive integer. 

The humps on the PDFs are most easily explained as 
an effect of grain "hopping" due to defects in the hexag- 
onal crystal structure. The humps can be observed only 
in relatively stiff systems where the hopping transport is 
not masked by the fluid vortical motion. This is clearly 
illustrated if we compare results for the stiff cluster (Fig. 
[5]) with those for the softer cluster. As seen from Fig. 3 
of Ref. in the latter experiment, PDFs are well de- 
scribed by Eq. d3J with p = 1.3 on the small time scales, 
and a transition to Gaussian shape is observed already 
at the scales of £ r ~ A. On very large scales, ~ 10A, 
new humps emerge (Fig. 4(a) of Ref. [lj|) which are due 
to large-scale vortical motion [3, G3] • 

In other experiments on quasi-2D laboratory dusty 
plasma similar results have been found. In studies of 
cold quasi-2D dusty-plasma liquids, Refs. @, & [HL [HJ 
reported non-Gaussian PDFs on scales below A. More- 
over, we find that the central portion of the PDF pre- 
sented in Fig. 3 of Ref. [ll[ can be fitted by Eq. (g]) with 
p = 1.15. On longer time scales, when particle diffuse to 
distances ~ A, the PDFs in those experiments approach 
a Gaussian shape. 

In the next subsection we show that the picture de- 
scried above is consistent with the results of simulations 
of the Langevin dynamics model. 



B. Langevin dynamics simulation results 

The effect of the hexagonal lattice structure is, as 
suggested above, most pronounced in the stiffer states. 
Langevin dynamics simulations with hard-wall confine- 
ment potential reveals this quite clearly. In Fig. G^d) 
one observes humps on the tail of the distribution up to 
~ 4A. Fig. [3^a), shows the PDF of horizontal displace- 
ments on time scales 5t, which is close to the collision 
time scale. On these time scales the grain motion is still 
ballistic, and the Gaussian PDF reflects the Maxwellian 
velocity distribution. On longer time scales the core of 
the distribution develops into a stretched Gaussian (see 
Fig. [3](c,d) and zoom-in (e,f)), and shoulders (and later - 
humps) develop on the tail. These non-Gaussian features 
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FIG. 3: (Color online) Langevin dynamics simulation of 
Yukawa system for the case of hard wall confinement. PDFs 
of horizontal displacements £r(r) = x(t + r) — x(i) (dots) in 
units of A, plotted for different time lags r as indicated in the 
upper corner. The full curve in (a) is a Gaussian fit and the 
full curves in the zoom-in figures (e,f) are fits by a stretched 
Gaussian given by Eq. Q with values of p as indicated in the 
upper left corner. 



are obviously a result of hopping of some grains from one 
potential minimum in the crystal lattice to another. 

It is interesting to compare these PDFs to those of the 
experiment displayed in Fig. [5J The magnified core of 
the distribution in Fig. GJd) corresponds roughly to Fig. 
[He). In the simulations, however, the statistics is better 
because of the larger number of grains in the cluster, so 
humps on the far tail can be observed in the zoomed-out 
image. These humps cannot be observed in Fig. HJc) due 
to limited statistics. They do appear, however at later 
times in the experiments, as shown in Fig. [He), when 
the width of the distribution exceeds A. This figure has 
not been reproduced in the simulations so far, because 
the simulations would have to run considerably longer. 

We have also analyzed the displacement PDFs for a 
sub-ensemble of particles in the core of the hard-wall 
cluster (sec Fig. [T{d)). These particles are completely 
trapped in their potential wells in the lattice. It turns 
out that in such a purely crystalline state the displace- 
ment PDFs on all time scales are Gaussian. 

Now let us discussed simulations performed for 
parabolic potential confinement. As discussed above, to 
simulate quasi-2D the motion in the vertical z compo- 




FIG. 4: (Color online) Langevin dynamics simulation of 
Yukawa system for the parabolic confinement case (stiff clus- 
ter). PDFs of vertical z-component of displacements £ 2 (t) = 
z(t + t) — z(t) (dots) in units of A, plotted for different time 
lags r as indicated in the upper right corner. The full curve 
in (a) is a Gaussian with the same standard deviation. 

nent is severely restricted. That results in the behav- 
ior of z-component of displacements similar to the hard- 
wall confinement case, as seen from Fig. |4j The PDFs of 
horizontal displacements, shown in Fig. [5l however still 
exhibit humps though much less pronounced because of 
the weaker confinement on (x, y) plane. As shown in Fig. 
[5jd) it approaches an overall Gaussian shape on longer 
time scales. 

In simulations of even softer states, though still with 
preservation of the regular lattice structure on the small- 
est time scales, stretched Gaussian distributions still ap- 
pear on the time scales for which the standard deviation 
of the displacements is less than A, but on longer time 
scales the PDF tends to develop humps on spatial scales 
much larger than A (Fig. |S](d)). Such features were also 
observed in the soft experiment (see [lij], Fig. 4(a)), and 
was attributed to the formation of vortical flows on spa- 
tial scales extending from A up to the cluster size. 

This soft state has also been run in a pure 2D configu- 
ration, i.e. without z-component, with results similar to 
those presented in Fig. [6] A preliminary fully 3D simula- 
tion of this state also yields qualitatively similar results, 
and this indicates that dimensionality is not a critical is- 
sue. The qualitative features described here are robust 
and universal for Langevin dynamics of Yukawa clusters 
with partial preservation of crystalline order. These stud- 
ies will be published in a forthcoming paper. 

Moreover, in these softer states, the system self- 
organizes into vortexes, or if stiffer, into rotating hexag- 
onal domains which continuously break up and merge. 
The effect, called "cooperative hopping", was first re- 
ported in 2D dusty plasma clusters in 0. Large-scale 
vortexes driven by stochastic forcing have been observed 
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FIG. 5: (Color online) Same as in Fig. [4] but now for the 
horizontal displacements 
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FIG. 6: (Color online) Same as in Fig. [5] but now for the soft 
cluster. 



experimentally [l4| and re pro duced by the code used in 
the present study in Ref. [131 . By plotting the particle 
positions and their pre- history over a suitable short time- 
interval can give an impression of the velocity-field of the 
vortical flow. If this is done in successive time frames 
one can construct movies of this velocity field, which are 
very instructive to watch. Such movies for examples of 
experimental and simulated flows can be found as auxil- 
iary material to Ref. [l7| • In Refs. [HI, [l4| one can find 
triangulation diagrams showing the occurrence of 5-fold 
and 7-fold defects associated with this hopping. 



IV. MEMORY EFFECTS AND 
INTERMITTENCY 



Information about the memory effects in the transport 
characteristics can be obtained from the variogram of 
particle displacements defined as 



S(t) = (N-r/St) 



N-r/St 

J'=i 



(5) 



where we have introduced the position increments 
£x ( r ) = x(j St + t) — x(j St) over the time lag r. For a 
fractional Brownian motion (fBm) the PDF of (r) is 
Gaussian, and S(t) = Dt 2H , where H is the Hurst ex- 
ponent. Normal diffusion corresponds to Gaussian PDF 
and H = 0.5 (Brownian motion). H different from 0.5 
is usually referred to as sub- or superdiffusion, and in- 
dicates violation of the central limit theorem, which im- 
plies that cither the displacements arc not identical and 
independent random variables and/or each of them has 
a heavy-tailed distribution with infinite variance. The 
latter is not the case in our experiments and simulations. 
The limit H — ► 1 corresponds to ballistic motion. 

In interpreting displacement variograms one has to re- 
alize that large-scale motions also influence the variogram 
on short time scales, so the Hurst exponent obtained on 
short-time scales does not necessarily reflect the nature 
of the diffusion of a particle relative to its nearest neigh- 
bors. Consider, for instance the extreme case of a frozen 
"flake" of dust grains which are completely fixed relative 
to each other, but that the flake itself is subject to trans- 
port characterized by a certain Hurst exponent. Then 
the variogram of every particle in the flake is charac- 
terized by this exponent, while their relative diffusion is 
negligible. In order to study this relative diffusion we 
would have to consider the evolution of relative position 
between neighboring particles. 

Information about H can also be obtained from the 
power spectral density of the particle displacement as a 
function of time. A power-law form of the spectrum, 
S(f) ~ f~ , corresponds to an fflm with H=(/3- l)/2. 
Thus normal diffusion requires (3 = 2. 

The main departure from Gaussianity observed in our 
displacement PDFs is a flatter (more heavy-tailed) ap- 
pearance, often well presented as a stretched Gaussian. 
A convenient measure of the flatness is the the kurtosis, 
which is the ratio between the fourth moment and the 
squared variance of the distribution. For the Gaussian 
the kurtosis is 3, and is higher for more heavy-tailed dis- 
tributions. In the theory of turbulence, the concept of 
intermittency is used for situations where the PDFs have 
higher kurtosis on small than on large scales (24j . In our 
context intermittency would require that the kurtosis of 
displacement PDFs decreases with increasing time scale. 
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A. Laboratory results 

For the laboratory data of the stiff cluster, the power 
spectral density, autocorrelation function and variogram 
of particle displacements are presented in Figs. 3(a,b) of 
Ref. [IH. In this experiment the PDF is non-Gaussian 
and self-similar on the short time scales, and the trans- 
port is superdiffusive with Hurst exponent H = 0.6 (Fig. 
6(a) of Ref. [Hj]). From results presented here, Fig. [21 it 
is apparent that on scales r 3> ta the displacement PDF 
approaches a Gaussian, so the kurtosis declines towards 3 
at longer time scales. This indicates some intermittency 
- the displacements are more bursty on small scales than 
on large scales. 

For the softer cluster, the power spectra are presented 
in Figs. 4(c,d) of Ref. [14| and the variogram in Fig. 3(b). 
The transport is superdiffusive with H = 0.84 on time 
scale r < ta and with H — 0.68 for r > ta- By watch- 
ing the movies of the particle motion one cannot avoid 
noticing that a common mode of motion the "minimal 
vortex" consisting of a stationary dust grain about which 
its six nearest neighbors rotate. This constitutes a nearly 
ballistic motion on the time scale ta, and might be the 
main contribution to the high value of H observed on 
scales t < ta- On longer time scales vortex motions on 
increasing scales are responsible for this memory-based 
supcrdiffusion. On the time scales one could follow the 
displacements in this experiment the PDF develops some 
large humps due to long-lived, large vortices, and the tails 
are cut off by the finite size of the cluster. Hence, this 
experiment cannot tell us about the fate of the PDF on 
very large scale in extended systems. 

In the experiments on cold dusty-plasma liquids re- 
ported in Ref. Q superdiffusive transport with H ~ 0.65 
is found for r < ta, and normal diffusion (H = 0.5) on 
larger scales. 



B. Langevin dynamics simulation results 

Variograms for different runs with varying stiffness are 
presented in Fig.fTJa). The variogram for the stiff core of 
hard-wall confinement is flat after a certain time (H — > 
0), and on average particles do not travel more than ~ 
0.1 A, and is not plotted in the figure. 

The curve represented by the • - symbols is represen- 
tative for particles from the looser boundary regions of 
hard- wall potential runs. Such particles are able to es- 
cape caging, but the transport is subdiffusive [H — > 0.3) 
even for long time scales. The curve has an "S-shape" 
typical for cases when an oscillatory component is present 
in the time-series. The same should be reflected in auto- 
correlation function of the signal. Indeed, in Fig. [8] we 
find a clear periodicity on the relevant time scale. 

For simulations of stiffer states with parabolic poten- 
tial (the curve represented by (A) - symbols) the caging 
of particles is still noticeable for r < 150 6t but for large 
times the Hurst exponent is higher, H — > 0.4. 




FIG. 7: (Color online) Variograms <t(t) = \/ S(t) in units of 
A, where S(t) is denned in Eq. ((5]l (top panel (a)). (•): 
Boundary region of hard-wall potential, corresponding to 
PDFs in Fig. [3] (a): Parabolic confinement the case of stiff 
cluster, corresponding to PDFs in Fig. [5] (T): Parabolic con- 
finement, the case of soft cluster, corresponding to PDFs in 
Fig. [6] Kurtosis as function of <t(t) (lower panel (b)). The 
symbols represent the same states as in panel (a). 



For simulations of soft states with parabolic potential 
the variogram (the curve with (T) - symbols the effect of 
caging is less prominent, but in the limit of large t we 
find also find H — ► 0.4 for this case. 

The general picture for all states simulated (except 
when caging is total) is that for t < ta the PDFs 
are stretched Gaussians (with humps superposed for the 
stiffer states), and the diffusion is slightly subdiffusive 
- more so for the stiffer states. These distributions ex- 
hibit heavier tails than Gaussians, and hence kurtosis 
(flatness) higher than 3. For t > ta the PDFs seem to 
tend more towards Gaussian again, at least when they 
are smoothed to remove the humps. This can be seen by 
plotting the kurtosis against the standard deviation as 
done in Fig. \T[b). It is observed that in general the flat- 
ness increases with time as long as the standard deviation 
< J { T ) = (C( T )) 1 ^ 2 ^ A, but have a clear descending trend 
towards the Gaussian value (3) as <j(t) exceeds A. For 
completely caged particles the kurtosis remains close to 
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dust grain in the form: 



T/Sf 



FIG. 8: (Color online) Displacement autocorrelation function 
(ACF) for the hard wall confinement case. The symbols rep- 
resent the same state as in Fig. [7J (a). The sampling interval 
used is 0.01 St. 



3, confirming that the PDF remains Gaussian for these 
particles. 

The concept of intermittency in a signal is loosely de- 
fined as the tendency of a signal to be more bursty on 
short scales than on long scales 24 1. This implies higher 
kurtosis on small than on large time scales. In our study 
this is the case on the "hydrodynamic" scales, r > ta, 
but not on the scales of caging, t<ta. 



THE ORIGIN OF NON-GAUSSIAN 
DISPLACEMENTS 



The Langevin thermostat used here models the interac- 
tion of the Yukawa many-body system with a heat bath of 
temperature T, and hence represents a canonical statis- 
tical ensemble. Thus the simulated system should relax 
to a thermodynamic equilibrium. With proper tuning 
of control parameters our system develops strong fluc- 
tuations, which in [l7j was shown to exhibit (imperfect) 
scaling. This makes it natural to think about this as 
equilibrium critical behavior. Since this behavior is ob- 
served only when the control parameters are in a range 
where the system is subject to some local positional or- 
der, we believe it is associated with a continuous solid- 
liquid phase transition. However, the small size and 
strong inhomogeneity of the systems simulated excludes 
direct comparison with theory [l[ and earlier simulations 
[25l ] for such transistions in 2D Yukawa systems. 

In extensive systems [(systems without long-range in- 
teractions)], thermodynamic equilibrium implies Gibbs- 
distributed microstates of the cluster. This, however, 
does not exclude non-Gaussian distribution of particle 
displacements on varying time scales. In order to obtain 
some insight into the origin of this non-Gaussianity it 
may be useful to write a Langevin equation for a single 



dr 

~dt 
dv 

It 



-[V<£(r) + V<f> st (r,t)} - v n v 



(6) 
dW , . 



where 4>{r) is the self-consistent (time-averaged) caging 
potential well that a partially caged dust grain experi- 
ences, 4> s t(r,t) is the additional fluctuating (stochastic) 
force field due to the Yukawa forces from the motion of 
the neighboring dust grains, v n is the dust-neutral colli- 
sion frequency, D = 2vT /m represents the diffusion coef- 
ficient due the stochastic force arising from collisions with 
neutrals, and W = (W x ,W y ,W z ) is the vector of three 
independent Wiener processes (Brownian motions). 

Let us first recall some results on free diffusion of a 
dust grain imbedded in a neutral gas, i.e. let us neglect 
the term in square brackets on the right hand side of 
Eq. ([7]). For a particle starting at rest at the origin at 
time t = 0, the equations can be integrated to yield [26| 
(r 2 (t)} — > (2T/mv n )t, for t > r„ = v~ x . On time scales 
t <C t„ the collision term v n w is unimportant, and the 
velocity is a Wiener process. On these time scales the 
position process is the integral of the Wiener process, 
and is ballistic. The variogram would yield H = 1. For 
t S> t„ the velocity increments will be independent and 
Gaussian (white noise), and the position increments (dis- 
placements) will be Gaussian and a Brownian motion, i.e. 
the slope is H = 1/2 in the variogram. 

If we introduce a stationary potential structure, but 
still neglect the stochastic force due to grain-grain col- 
lisions, there is a well-defined Fokker-Planck equation 
formulated in 6-dimensional phase space associated with 
Eq. (J7J, which features the Gibbs distribution P(r,v) ~ 
cxp[— (^mv 2 + Q<j)(r)) /T] as a stationary solution, and in 
configuration space we have the Boltzmann-distribution 
p(r) ee / P(r, v) d 3 v ~ exp [-Q0(r)/T]. 

If the caging potential barriers are much higher than 
the thermal energy of the grains, the cluster is in a com- 
pletely frozen crystalline state, there is no diffusive trans- 
port beyond the mean interparticle distance. 

The omission of the stochastic force due to grain-grain 
collisions is only valid on time scales r <C r c , where r c is 
essentially the time it takes for a grain to bounce off the 
walls in its cage. The typical situation in our experiments 
and simulations is that 5t ~ t c <C t„. On time scales 
t < t c <C r„, as explained above, the motion is ballis- 
tic. If the velocities are Maxwcllian due to the Langevin 
thermostat, the displacements £(r) ~ vt will be Gaus- 
sian distributed for r < t c . This is what we observe in 
the simulations. In the experiments discussed in this pa- 
per we observe non-Gaussian displacements already for 
t = 8t, but this is because the sampling time in the 
experiment is somewhat longer than the grain-grain col- 
lision time. This is supported by the results reported in 
[l8| , where Gaussian displacements were found for r = 5t 
in an experiment where particle positions were sampled 
at 500 Hz. 
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On time scales r > t c a single-particle stochastic equa- 
tion makes sense only if the force due to grain-grain inter- 
action is represented as a stochastic source term. If the 
issue under consideration is grain transport, the stochas- 
tic equation can be simplified by averaging the equations 
over a time scale of the order r n , and obtain an equa- 
tion for the averaged position. This justifies to neglect 
the inertial term dv / dt in Eq. ([7]) (the Smoluchowski ap- 
proximation), which then reduces to the 



dr 

dt 



Q 



V0(r) 



(8) 



where F st represents the stochastic force due to collisions 
with neutrals as well as other dust grains. A problem is of 
course that the interaction F st contains long-range spa- 
tiotemporal dependencies and the statistical distribution 
of momentum transfers may be non-Gaussian, and hence 
cannot reliably be represented by a white noise term. 
One could certainly make assumptions about the statisti- 
cal nature of the interaction, but the standard derivations 
of the associated Fokker-Planck equation are no longer 
valid since the position is no longer a Markov process. 

Another approach would be to assume that the time- 
averaged potential <j) is constant in the interior of the 
cluster, and to incorporate the the effects of the dust-dust 
interaction in the friction term and the stochastic term 
in Eq. ([?]). The resulting stochastic equation would take 
the form (for simplicity we write the one-dimensional ver- 
sion): 



dv 
~dt 



K{v) + y/D(v) 



dW 
~dT' 



(9) 



where the drift term K(v) represents a generalized fric- 
tion due to combined action of the dust-dust collisions 
and the dust-neutral collisions, and the second term on 
the right represents the stochastic part of these interac- 
tions. In the following we shall assume that the latter 
contains no long-range memory and can be represented 
as Gaussian white noise. An equation of this form was 
studied by Kaniadakis and Quaraty (2?| , and by Borland 
[28| in the context of determining the conditions K(v) 
and D(v) must satisfy for the resulting stationary veloc- 
ity PDF to reduce to the Tsallis distribution^. If the 
stochastic e quat ion is interpreted in the sense of ltd (see 
for instance [2y]), the associated Fokker-Planck equation 
for the velocity distribution P(v,t) is 



dP(v, t) 
dt 



^\K{x)P{v,t)\ + ~[D(v)P{v,t)]. 



dv 



(10) 

If the drift term has a form that confines the particles 
to a range with finite velocities (which will be the case 
for a physically meaningful friction term) the stationary 
solution to Eq. (fit))) can be written as 



P(,)=iVexp(/^[A»-i 



1 dD(u) 



On 



du 



(11) 



where N is a normalization constant. In the case of 
linear friction K(v) = —v c v and constant diffusion co- 
efficient D(v) = D Eqs. © and (fTO")) represents the 
well-known Ornstein-Uhlenbeck problem, and the solu- 
tion (fTTj) is the resulting Maxwell velocity distribution. 
Borland [28[ shows that a general criterion for the solu- 
tion (fTTj) to reduce to the Tsallis distribution [29[ 

P g (v)=Z- 1 [l-/3(l-q)v 2 ] 1 ^-^ (12) 
is that {K(v), D(v)} satisfy the relation 



D(v) 



K{v) 



ld_ 

2di 



D(v) 



0v 



(13) 



Here (3 — \/2kT is the inverse temperature and q is the 
non-extensivity parameter. For q = 1 the distribution 
reduces to the Maxwellian, and the family of functions 
that gives this conventional extensive Gibbs-Boltzmann 
statistics is given by the relation 



K{v) 



1 dD 
2~d~u 



(3Dv. 



(14) 



For constant diffusion coefficient this reduces to K(v) = 
— v c v, where v c satisfies the relation v c = D/2kT, which 
is a version of the fluctuation-dissipation theorem. 

Borland emphasizes that Eq. (fl3|) implies that there is 
a whole family of functions {K(v), D(y)} which give rise 
to Tsallis statistics, but it should maybe also be men- 
tioned that there exists an infinitely larger family of func- 
tions that do not. Thus, the ubiquity of Tsallis statistics 
in this context relies on the existence of ubiquitous phys- 
ical situations where Eq. p3[) is fulfilled. The analysis of 
our experiments and simulations shows that the displace- 
ment PDFs are well described by stretched Gaussians 
with p w 1.5 on time scales in the range r c < r < ta, 
and by Gaussians in the range r > ta . In Fig. [5] we 
demonstrate that the stretched Gaussian with p = 1.5 
is practically indistinguishable from the Tsallis distribu- 
tion with q = 1.1, so there is no information in the ex- 
perimental data to prefer one to the other. On the other 
hand the stretched Gaussian provides a very simple and 
physically interpretable generalization of the Ornstein- 
Uhlenbeck problem, as will be described in the following. 

Let us consider the Ito-Langcvin equation Eq. ([9]) with 
constant diffusion coefficient D, and generalized friction 
term 



K(v) = -sg\i(v)v c \v\ 



(15) 



The physical rationale for such a sub-linear friction term 
is simply to account for the reduced effective dust-dust 
friction introduced by those fast particles undergoing co- 
operative hopping. It has been shown 0, [|| that this 
hopping takes place along string- like structures of defects, 
and hence it is more natural to model it as a modification 
of the directed friction force than as modification of the 
diffusion coefficient. 
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FIG. 9: (Color online) The red curve is a log-plot of the 
stretched-Gaussian function exp(— Bx p ) vs. x 2 for B — 1.2 
and p = 1.5. The blue curve is a corresponding plot of the 
Tsallis distribution [1 + (q - l)x 2 ]- (q - ir 1 for q = 1.1. The 
values for p and q are typical for the experimental results 
discussed in this paper. 



In Eq. [15] the case p = 2 corresponds to the standard 
Ornstein-Uhlcnbcck problem, and p = 3/2 would corre- 
spond to a friction that grows as the square root of the 
velocity, rather than proportional to the velocity. The 
stationary solution to the associated Fokker-Planck equa- 
tion given by Eq. (jlip reduces to the stretched Gaussian 



P p (v) = N p exp[-{2 Vc /p)\v\n]. 



(16) 



At this point we should keep in mind that the stretched 
Gaussian Eq. (fj| found in our data are PDFs for the 
position displacements at varying time lags r, not the 
velocity PDFs. In order to find the displacement PDFs 
we have integrated numerically the stochastic equation 
with the generalized friction term (and p = 1.5) for 
an ensemble of initial velocites drawn from the strctched- 
cxponential distribution Eq. (|16p . Then we have found 
the position trajectories by integrating v(t), and formed 
the displacement PDFs for different time lags. The rather 
obvious result is that the displacement PDF remains 
close to the stretched Gaussian for r ~ r c = v~ x and 
converges to a Gaussian for r ^> r c . The reason for this is 
that the displacements are near ballistic for r < r c , and 
hence reflect the velocity distribution. On time scales 
greater than r c the velocity process is de-correlated and 
stationary, and successive displacements can be consid- 
ered as independent, identical random variables. Dis- 
placements £ r , for t ^> t c , can therefore be considered 
as a sum of identical, independent random variables, and 
thus the conditions of the central-limit theorem are satis- 
fied. Thus the displacement PDF should be Gaussian and 
the diffusion should be normal (Brownian motion). The 
variograms produced on the basis of the computed PDFs 
shows (obviously) that H = 1 for r <C t c and H = 0.5 for 
t ^> r c , but also that there is a large transitional regime 
of time scales r c < r < 50r c where the slope corresponds 
to Hurst exponents in the range 0.5 < H < 1. 

This extremely simplified stochastic model reproduces 



the stretched exponentials observed for r c < r < ta , and 
the convergence to Gaussian for t ta observed in our 
simulations. Since the sampling time in the simulations 
are of the order of t c the variograms in Fig. [7] do not 
show the ballistic regime where H = 1 , but the behavior 
of the variograms in the regime t ta can be inter- 
preted as a slow convergence towards a Brownian motion 
(H = 0.5). It is also noteworthy that the simulations 
in [H, Q gave rise to similar variograms when the system 
was in the hexatic phase between a solid and a liquid. 
Those simulations were performed with periodic bound- 
ary conditions, so the effects of the finite systems size 
and the system inhomogeneity did not influence the vari- 
ogram at large times as they do in our simulations. Since 
the simulation in Q also displays a leptokurtic distribu- 
tion on the scales r < r c , and a transition to a Gaussian 
for t > t c , we conclude that the stochastic model, and 
the associated Fokker-Planck equation, yields a reason- 
able lowest-order description of anomalous diffusion in 
2D many-body systems in the hexatic phase. 

The linear Fokker-Planck approach discussed above 
(including the cases discussed by Borland) implies Gaus- 
sian displacement distribution in the limit r — > oo. Hence 
it cannot be used to describe anomalous diffusion in 
the long-time limit. Another formulation that allows 
non-Gaussian displacement distributions and anoma- 
lous transport is the nonlinear Fokker-Planck equation 

tun, 



dP(x, t) 
di 



(17) 

In order to describe anomalous displacement diffusion x 
should now be interpreted as a spatial variable, and we 
should focus on time-dependent solutions. For this par- 
ticular nonlinear diffusion term, and a linearly varying 
drift term K (x) — k\ + k^x^ such a time-dependent so- 
lution takes the form of a Tsallis distribution, Eq. (|12[) . 
where (3(t) and Z q {t) now are time dependent. In the 
case of free diffusion k\ = ki = 0, the variance l/(3(t) of 
the distribution has the time-asymptotic dependence 



t a , a = 



(18) 



provided q < 5/3. For q > 1 the distribution has a power- 
law tail P q (x) ~ x 2 /' 9-1 ), and for q > 5/3 the tail is so 
heavy that the variance does not exist. For 1 < q < 5/3 
we have power-law tail and a > 1, implying superdiffu- 
sion. For q = 1 the distribution reduces to a Gaussian 
and we have normal diffusion, and for q < 1 the Tsal- 
lis distribution has a cut-off in the tail (P q (x) = for 
x > [1 — /3(t)(l — g)]~ 1//2 ]), and the transport is subdif- 
fusive. Thus the Tsallis relation (fT5|) implies that distri- 
butions with tails heavier than the Gaussian (leptokur- 
tic) must be associated with superdiffusion. Our sim- 
ulations are subdiffusive for r > r c , but the PDFs are 
leptokurtic and well described by the Tsallis distribution 
with q w 1.1, which according to Eq. (fT5|) should give a 
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weakly supcrdiffusivc transport. This inconsistency indi- 
cates that the nonlinear Fokker-Planck equation is not a 
good description of the transport observed in our simu- 
lations. 

There are also considerable difficulties associated with 
the physical interpretation of Eq. (fT7|) for our case. Bor- 
land [3l| attempts to provide such an interpretation by 
considering D oc P x ~i as a P-dcpcndent diffusion coeffi- 
cient, that could replace D in the stochastic equation © 
(since we are considering spatial diffusion, the variable v 
in that equation should now be interpreted as a spatial 
variable, so we are really dealing with the Smoluchowski 
equation ((8|)). If P(x,t) could be interpreted as a den- 
sity, there would be no problem to envisage a density- 
dependent diffusion coefficient, but the density of dust 
grains in our system is almost constant in space and time, 
so a spatial diffusion equation for density does not make 
sense. P(x, t) will have to be interpreted as a probablity 
density for an ensemble of independent particle trajecto- 
ries starting out at positions selected from a prescribed 
initial distribution P(x,0) in different independent real- 
izations of the experiment. Since we are interested in 
the evolution of particle displacements a natural choice 
is to consider an ensemble of particles starting from the 
same postition xq, i.e. P(x, 0) = 6(x — xq)- The quantity 
P(x,t) obviously depends on the choice of initial pos- 
tion xq, and therefore D oc P 1 ^ q (x,t) cannot be a local 
physical property of the medium at the point at (x,t). 
It is rather a quantity that is valid only for the parti- 
cles belonging to this specific ensemble (i.e. the particles 
starting out at the position Xq at t = 0), and hence such 
a model would imply that the diffusion properties of a 
particle depend very strongly on its prehistory, i.e. there 
must be strong long-range memory in the system. 



VI. CONCLUSIONS 

In this paper we have studied anomalous transport 
properties; non-Gaussian particle displacement statistics 
and memory effects, in a Langevin-dynamics model, and 
we have contrasted these simulations with experimental 
results from quasi-2D dust plasma clusters in soft states. 
We have demonstrated that anomalous transport is ubiq- 
uitous on time scales t c < t < t\ provided the system 
is not residing in the pure crystalline or liquid phases. 
Moreover, we have shown that this complex dynamics 
emerges from stochastic forcing and dissipation provided 
by the collisions with the neutral gas, and that no other 
plasma effects than the static Debye screening arc re- 
quired. 

From the experiments and the Langevin dynamics sim- 
ulations we observe that Gaussian displacement distri- 
butions are restricted to the following cases: (i) on time 
scales less than the collision time r c in all systems, (ii) 
on all time scales for systems that are so stiff that no 
hopping takes place (pure crystal) , (iii) on all time scales 
for systems where no caging occurs (pure liquid or gas), 



(iiii) for time scales longer than ta in soft systems. 

A result that seems to be robust for both stiff and soft 
states, is that PDFs evolve from a leptokurtic (heavy- 
tailed) shape for small time-lags t<ta towards a Gaus- 
sian shape for larger r. This feature is usually interpreted 
as a signature of intermittency - a spikyness in the flow- 
field which gives rise to more leptokurtic PDFs of the dis- 
placements on the small scales. Note that intermittency 
does not prevail to spatial scales less than A. A simple 
way to model this intermittency feature is via the singlc- 
partice stochastic equation @ , and the associated linear 
Fokker-Planck equation (JTUJ) , with a non-linear drift term 
given by Eg. (| 151) . 

In Ref. [141 ] it is pointed out that the dust mono- 
layer should be perceived as a viscoelastic medium, where 
the elastic properties dominate on the smaller scales. 
Thus the transition from stretched-Gaussian to Gaussian 
PDFs could correspond to the transition from elastic to 
viscous flow. One future goal would be to perform exper- 
iments and simulations on sufficiently large spatial and 
temporal scales to characterize the transport associated 
with homogeneous, turbulent flow in such systems. With 
such extensive simulations one could clarify, for instance, 
if energy will cascade to vortices on infinitely large scales 
in an infinite system, or if there is an upper limit to the 
vortex size related to the characteristic length over which 
elasticity dominates over viscosity. 

The rather imprecise qualifier "soft" has been em- 
ployed in this paper to characterize 2D dusty plasma 
clusters which preserve local hexagonal order, but ex- 
hibit loss of such order on larger scales. Materials dis- 
playing such properties are also classified as "hexatic" , 
and the continuous transition from the crystalline to the 
liquid state might be considered as a continuous phase 
transition In (l7j we perform a particular avalanche 
analysis of the experimental soft cluster and a Langevin- 
dymamics simulation. In this analysis we define con- 
nected space-time regions where the kinetic energy of the 
grains exceeds a prescribed threshold, and define the ac- 
tivity in these regions as avalanches. The statistics of 
avalanche sizes and durations tend to obey power laws 
with similar characteristic exponents for experiment and 
simulation, which are properties characteristic for crit- 
ical phenomena and higher order phase transitions 0]. 
Such phenomena have been studied in strongly idealized 
models such as the Ising model [§, 33] and the XY model 
[HI , but the Yukawa-Langevin model studied in this pa- 
per exhibits, in spite of its simplicity with respect to the 
physical interactions, the prospect of quantitative pre- 
diction of dynamics in real physical systems where the 
microscopic state of the system can be tracked to the 
smallest detail of physical interest. 
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